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Abstract 

Compressible stability of growing boundary layers is studied by numerically 
solving the partial differential equations under a parabolizing approximation. The 
resulting parabolized stability equations (PSE) account for non-parallel as well as 
nonlinear effects. Evolution of disturbances in compressible flat-plate boundary 
layers are studied for freestream Mach n um bers ranging from 0 to 4.5. Results indi- 
cate that the effect of boundary-layer growth is important for linear disturbances. 
Nonlinear calculations are performed for various Mach numbers. Two-dimensional 
nonlinear results using the PSE approach agree very well with those from direct nu- 
merical simulations using the full Navier- Stokes equations while the required com- 
putational time is less by an order of magnitude. Spatial simulations using PSE 
have been carried out for both the fundamental and subhaxmonic type breakdown 
for a Mach 1.6 boundary layer. The promising results obtained in this study show 
that the PSE method is a powerful tool for studying boundary-layer instabilities 
and for predicting transition over a wide range of Mach numbers. 
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I. Introduction 

The subject of compressible boundary-layer stability has attracted a great deal 
of interest in the past few years due to its importance in understanding the on- 
set of transition in high-speed flows and providing some theoretical background 
for laminar flow control (LFC) techniques (Malik, 1990a). Most investigations of 
compressible linear stability (e.g., Mack, 1969, 1984) have employed what is known 
as the “quasi-parallel” approach whereby the growth of the boundary layer is ig- 
nored and the linearized Navier-Stokes equations are reduced to ordinary differential 
equations (ODE) by assuming a wave-like disturbance of the form 

#*,»,*,<) ~»(y)e i( " + * , -" <) (1) 

where x, y, and z are the streamwise, wall-normal, and spanwise coordinates, respec- 
tively; a and j3 are the corresponding wave numbers, u> is the disturbance frequency 
and 'it represents the disturbance shape function. The linear ODE’s along with the 
homogeneous boundary conditions constitute an eigenvalue problem of the form 

a = a(u, P ) ( 2 ) 

which can be solved by standard eigenvalue techniques. The imaginary part of 
a gives the disturbance growth rate and a small disturbance is expected to grow 
provided < 0. For a given flow, this eigenvalue approach can be applied “locally” 
at various locations along the body in order to obtain an idea about overall growth 
of disturbances and to correlate with transition location using empirical methods 
such as the e N method. 

The effect of non-parallel flow on boundary-layer instability has been studied 
by Gaster (1974), Saric and Nayfeh (1975), Gaponov (1981), and El-Hady (1991). 
In the multiple-scales method used by the latter three authors, the disturbances 
are decomposed into a slowly- varying shape function and a rapidly-oscillating wave 
part. Both parts are represented as functions of a fast-scale variable (x) and a slow- 
scale variable (x = ex, with e = 1 /R). With these assumptions the governing PDE’s 
are reduced to a set of ODE’s by neglecting terms of order equal to or higher than 
e 2 . In conjunction with the solvability condition, the analysis yields non-parallel 
corrections to the eigenvalues computed by the quasi-parallel theory. Just like the 
traditional linear theory, the multiple- scales approach can only be applied locally 
for a given problem. 

Apart from the “local” methods described above, the evolution of disturbances 
in a given flowfield may also be computed numerically by solving the governing 
partial differential equations (PDE’s) without resorting to the eigenvalue approach. 
The effect of boundary- layer growth and other history effects associated with ini- 
tial conditions and variation in wall temperature, for instance, can be properly 
accounted for. This was done for the Gortler vortex problem by Hall (1983) and 
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Spall and Malik (1989). Denier et al. (1991) solved the “receptivity” problem to 
provide the inflow conditions for the PDE’s and were able to show how Gortler 
vortex structure develops from a discrete roughness site. 

The governing PDE’s for the Gortler problem are parabolic and thus the solu- 
tion can be obtained by direct marching provided the initial conditions are known. 
However, the governing equations for Tollmien-Schlichting (TS) and inviscid type 
disturbances are elliptic and their solution cannot be obtained by simple march- 
ing methods. In addition, the numerical solution of these PDE’s requires proper 
outflow boundary conditions which is a nontrivial task. However, we note that for 
boundary-layer type flows which are of interest here, the equation set is only weakly 
elliptic along the dominant flow direction. Therefore, with appropriate simplifica- 
tions, one could “parabolize” these stability equations and avoid the difficulties 
associated with the downstream boundary conditions. 

From a physical view point, the streamwise ellipticity arises from the upstream 
propagation of acoustic waves and the streamwise viscous diffusion. To render 
the stability equations parabolic, one must devise a way to suppress, but without 
compromising the essential physics, this upstream propagation. One way to derive 
the parabolized stability equations (PSE) is to borrow ideas from the multiple- scales 
approach and decompose the disturbance into a rapidly- varying wave- like part and a 
slowly- varying shape function. The ellipticity is retained for the wave part while the 
parabolization is applied to the shape function. The resulting PSE can be solved by 
marching along the streamwise direction. The technique can be used to study both 
the linear and nonlinear evolution of convective disturbances in growing boundary 
layers. Global or absolute instabilities can not be studied by this approach. This 
parabolizing procedure has been used recently by Bertolotti et al.(1992) for Blasius 
flow. 

The objective of this research is to study compressible boundary-layer stability 
and transition. We employ parabolized stability equations for linear and nonlin- 
ear development of disturbances in a compressible boundary layer. The nonlinear 
calculations are carried all the way to the transition stage for supersonic flows. In 
section II, we formulate the problem while the numerical procedure used to solve 
the governing equations are given in section III. The results and conclusions are 
presented in section IV and V, respectively. 


II. Problem Formulation 

The evolution of disturbances in compressible boundary layers is governed by 
the compressible Navier-Stokes equations 

| + v-(^) = o 
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( 3 ) 


P [^~ + (Y ■ V)?] = - V P + V[A(V • V)] + V • [p(VV + VI? r )] 
at 

pcA^ + (V ■ V)T] = V . (*VT) +^ + (V.V)p+<i 

where V is the velocity vector, p the density, p the pressure, T the temperature, c p 
the specific heat, k the thermal conductivity, p the first coefficient of viscosity, and 
A the second coefficient of viscosity. The viscous dissipation function is given as 

$ = A(V ■ Vf + |[W + W T ] 2 . 

£ 

The equation of state is given by the perfect gas relation 


p — p^RT 

and the steady state solution of the basic flow can be derived by invoking the 
boundary-layer assumption. 

In this research, we formulate the compressible stability problem in Cartesian 
coordinates for the flat-plate geometry, although the theory itself can be easily 
extended to axisymmetric bodies and infinite swept-wing flows. The Cartesian 
coordinates are denoted by x, y, and z to represent the streamwise, wall-normal, 
and spanwise directions, respectively. All the lengths are scaled by a reference length 
l, velocity by u e , density by p e , pressure by p e a 2 , time by l/u e , and other variables 
by the corresponding boundary-layer edge values. The basic flow is perturbed by 
fluctuations in the flow, i.e. the total field can be decomposed into a mean value 
(boundary-layer solution) and a perturbation quantity 

u = u + u' , V = v + v' , w = w + w' 

p-p + p', p = p + p', T = T + T' (4) 

p — p -f* p! , A = A + A', k — k T k . 


Substituting Eq.(4) into the Navier-Stokes equations given by Eq. (3) and sub- 
tracting from the governing equations corresponding to the steady mean flow, and 
using the equation of state, we obtain the governing equations for the disturbances 
as 


dd> d6 dd> dd> 

+ A + B-£ + c^- + D<t> 

at ox ay az 


. V &± + v ** 

- y t.x * n i y xy 


dx 2 


dxdy 


d 2 d> 

_i_ v Z 
+ V V9 n.,2 


+ v; 


d^_ 

dxdz 


+ V t 


y* 


d 2 <t> 

dydz 


dy 2 

+ V t 


d 2 cj> 

dz 2 


( 5 ) 


where <j) contains the disturbance vector and is defined as 


$ - (p’,u',v',w',T') 


*\T 
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Matrices T, A, B, C, D , V**, K y , V yy , V**, V yr , and V** are Jacobians of the 
corresponding total flux vectors and are composed of a linear part with only mean 
flow quantities (denoted by superscripts l) and a nonlinear part which contains 
perturbation quantities (denoted by superscripts rt): T = T l + T n , A = A 1 + A n , 
etc. We note here that matrices T, A, R, C, D have contributions from both inviscid 
and viscous terms, and thus contain terms of order one and of order I/Rq (Rq is 
the reference Reynolds number Rq = u e l/u e ); while matrices V xx , V xx , V xy , V yy , 
V xz , V yr , and V zz are solely due to viscous diffusion and are of order I/Rq. 

To facilitate our discussion on the relation between linear and nonlinear dis- 
turbances, we rearrange Eq. (5) in the following form 




xx dx 2 Xy dx dy 


- V, 


i d 2 4> 

dy 2 


yy 


- VI 


dxdz 


7 l 

yz 


d 2 4> 

dydz 


V> = F » 

zz dz 2 


( 6 ) 

where the left hand side contains only linear operators operating on the disturbance 
vector and the right-hand-side forcing vector F n is due to nonlinear interaction and 
includes all nonlinear terms associated with the disturbances. The right hand side 
is given as 

F n = — r n ^ - A n ^- - B n ^- - C n ~ 
at dx dy dz 


- D n 4> + V n — + V n + V n ^ 

> r X T ^ o 1 r rw a r\ i r j r 


dx 2 


xy 


dxdy 


yy 


dy 2 


( 7 ) 


+ V" + v» + K" 

xz dxdz ^ yz dydz zz dz 2 ' 


In the incompressible limit, jF” contains quadratic nonlinearities; while, for com- 
pressible flows, cubic and higher-order nonlinearities are present. For small distur- 
bances, F n can be neglected and thus Eq. (6) reduces to the linearized Navier-Stokes 
equations 

dy dz 


T l<tt + a 1 — 

dt + dx 


+ B l 


-Vi 


d 2 <f> 


xx dx 2 

- V ‘£±. 

yz dydz 


-V 


- V. 


d 2 <f> 

xy dxdy 

d 2 4> 


dz 2 


- V. 


= 0. 


d 2 <j> 

dy 2 


-Vi 


ay 

dxdz 


( 8 ) 


The governing PDE’s of the disturbances, Equation (6), is hyperbolic in time 
for the convection terms (inviscid part). When we consider only the spatial deriva- 
tives, Equation (6) is elliptic in the streamwise direction due to two reasons. First, 
the streamwise viscous term V xx allows any disturbances to be diffused upstream. 
Second, and more importantly, the convection term in the streamwise direction 
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makes the upstream propagation of acoustic waves possible. The latter can be bet- 
ter understood by considering the linearized version of the inviscid equations and 
using the method of characteristics (MOC) theory. Since the inviscid part of Eq. 
( 8 ) is hyperbolic in time, the corresponding slopes of characteristic lines in the x — t 
plane (which determine the direction of propagation) can be found by solving the 
following eigenvalue equation 


| A 1 - A c r'| = 0. 

Negative eigenvalues imply the wave is propagating from downstream to upstream 
and vice versa. The eigenvalues of the above equation are 


A c = u, u, u, u + c, u — c 


where c is the speed of sound. For boundary-layer flows of interest in this study, the 
first four eigenvalues are always positive, while the last eigenvalue ( u — c ) can be 
either negative or positive depending upon the local Mach number ( M x = ujc). For 
subsonic flows, this quantity is negative throughout the whole flowfield, therefore, 
the equation set ( 8 ), and thus ( 6 ), is elliptic. For supersonic flows, the ellipticity 
only arises inside the subsonic layer adjacent to the wall. 

Based upon the above discussion, one way to “parabolize” the PDE’s given by 
Eq. ( 6 ) and make the marching solution feasible is to neglect the viscous diffusion 
terms along the streamwise direction and prohibit the upstream wave propagation 
either by dropping the left-running characteristics (associated with the eigenvalue 
u-c ) (Chang and Merkle, 1989) or suppressing some part of the streamwise pressure 
gradient, as it is done in the Parabolized Navier-Stokes (PNS) approach (Vigneron 
et al., 1978). For the stability equations, the upstream wave propagation can be 
suppressed by either dropping the characteristic equation associated with the eigen- 
value u-cor multiplying the streamwise pressure disturbance gradient dp' /dx by 
a parameter fl given by 


7 Ml 

i+( 7 -i wr 

= 1, M x > 1 


M x < 1 


( 9 ) 


where 7 is the ratio of specific heats. These parabolizing procedures are quite 
effective for the PNS approach and yield solutions which compare favorably with 
those obtained by the full Navier-Stokes equations provided a large portion of the 
flow is supersonic and only steady state solutions are of interest (Vigneron et al., 
1978; Rubin, 1981). The advantage, of course, is the significant reduction in the 
computational cost due to the marching solution. 

For compressible stability problems, the disturbances are essentially unsteady 
waves propagating across the whole boundary layer and the amplitudes of these 


5 



waves reach their maxima near the critical layer located between the wall and the 
boundary-layer edge. These instability waves undergo a “fast oscillation” (phase 
change) as they evolve along the flow direction. Direct application of the parab- 
olizing procedure used in the PNS approach for mean flow computations to our 
governing stability equations would not capture the flow physics due to the sup- 
pression of the wave propagation along the left-running characteristics. Therefore, 
an alternative procedure must be devised. 

Linear PSE 


As mentioned previously, one way to “parabolize” the governing PDE’s is to 
first decompose the disturbances into a fast-oscillatory wave part and a slowly- 
varying shape function. We keep the ellipticity for the wave part while parabolizing 
the governing equation for the shape function. Following the lead of the non-parallel 
linear stability theory, we assume that the disturbance vector <f> for an instability 
wave with a frequency u and a spanwise wave number ft (assume the wave is periodic 
in both the temporal and spanwise directions) can be expressed as 

<f>(x,y,z,t) = ^(x,y)e'^’’o Q( -^ di+ ^ z (10) 


where x is the fast-scale variable, a(x) is the corresponding streamwise wave number 
and $ is the “shape function” vector given by 

= (p,u,v,w,T) t . ( 11 ) 


As compared to Eq. (1), the shape function 'F is now a function of both x and y 
due to the growth of the boundary layer and the wave number a is a function of x 
to account for the growing boundary layer. For simplicity, we now restrict ourselves 
to the linear case, i.e., only a single disturbance mode (a;, /3) is considered and the 
nonlinear effect will be included later on. Substituting Eq. (10) into the linear 
stability equation (8), we have the following equation for the shape function 


DV + A^+B^- 

ox ay 


dx 2 


uJLuy 


( 12 ) 


where the vectors 


D, A and B are defined by 


D = -iuT 1 + D l + iaA 1 -|- ipC 1 

- (<^ - aW, + affVi, + pv'„ 
A = A 1 - 2 iaVj, - i/3V‘ z 
B = B'-iaVl,-i0V‘,. 
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In the quasi-parallel linear theory where “normal-mode” analysis is employed, 
the shape function 'F is assumed to be a function of y only (city /dx = 0); therefore, 
Equation (12) reduces to the following system of ODE’s 


Lgty = 0 


(13) 


where the operator Lq is given by 


L o 


= D + 6 T y- V » 


dy 2 


and the elements , of matrices D, B and V yy axe evaluated by assuming parallel 
mean flows (u = 0 and da/dx = 0). The above ODE’s in conjunction with homoge- 
neous boundary conditions then constitute an eigenvalue problem described by the 
dispersion relation given in Eq. (2). 

Unlike the normal-mode analysis described above, the decomposition (10) ex- 
hibits some extent of non-uniqueness between the distribution of the wave part and 
the shape function part. In the PSE approach, we choose a complex wave num- 
ber a and construct a decomposition such that the change of shape function 'F 
along the streamwise direction x is of order 1 / Rq and the second derivative of '3/ 
(d 2 ty / dx 2 ) is negligible. With this assumption and after neglecting all terms of 
0(1/ Rq), Equation (12) reduces to 


~dty 
Dty + A — 

ox 


- 3 * 

dy yy dy 2 


(14) 


Equation (14) describes the evolution of the shape function 'F and is “nearly” 
parabolic in the sense that second derivatives in x are absent and the elliptic effect 
associated with the wave part is absorbed in matrices D , A and B. For instance, 
the disturbance pressure gradient dp 1 [dx, which is responsible for the upstream 
influence, can be written as 


dp/ 

dx 


„ dp . i( f a{x)dx+j3z—uit) 

— (tap + 77— )e 
dx 


The contribution of the wave part (iap) is absorbed in the source term D’F and does 
not contribute to the upstream influence of the governing equations of the shape 
functions, Eq. (14). However, the pressure gradient shape function dpjdx associ- 
ated with the left-running characteristic (for subsonic flows only) is still present in 
the x derivative term. The existence of this term allows upstream influence in Eq. 
(14). For stationary Gortler vortex problem; a = 0 and dpjdx drops out, Eq. (14) 
reduces to the parabolic equations solved by Spall and Malik (1989). 

For supersonic boundary layers, a large portion portion of the flow possesses 
only downstream characteristics, our numerical results have shown that with a 
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properly chosen value of a (see discussion below) most of the upstream influence 
is accounted for in iap and the elliptic effect associated with the pressure gradient 
shape function, dp/dx, is negligible. To make Eq. (14) truly parabolic and enable 
a stable marching procedure for subsonic flows, we multiply dp/dx by a constant 
Q defined in (9) (in the incompressible limit, this is equivalent to setting dp/dx to 
zero). While, this is formally true only in special cases (e.g. Gortler vortex problem), 
the approximation yields solutions which compare very well with accurate results 
from full Navier-Stokes equations (Joslin et al., 1992). This is because most of the 
ellipticity is captured in the iap term. 

For incompressible flows, one can use the vorticity-streamfunction formulation 
for two-dimensional flows or use other formulations derived by eliminating the pres- 
sure from the momentum equations, as is done by Bertolotti et al. (1992). In these 
approaches, neglecting second and higher streamwise derivatives of the dependent 
variables inherently suppresses some part of the streamwise pressure gradient, and 
consequently prohibits the upstream propagation of information. 

We now describe the strategy to update the streamwise wave number in or- 
der to make the marching scheme well-posed. The evolution of shape functions is 
monitored during the process of marching and a is updated by local iterations at a 
given x according to the change in The updating procedure is described herein. 
At a given location x \ , we assume that the streamwise wave number is given by Oj 
and the total disturbance in the vicinity of Xi can be expressed as 


<f>(x,y,z,t) = $(x,y)e 


[( r a\dx+ fiz— u><) 

j XI 


(15) 


The change of the shape function 4' can be approximated by the following Taylor 
series expansion truncated to the first order 


$(x,y) = + -^-(x - xi) ... 

where 'kj is the shape function at x — X\. To an accuracy of 0(x — xx), the above 
equation can be further expressed as 


p i <i*i fT . 

$(x,y) = ^ . (16) 

Substituting (16) into (15), we have the “effective” wave number in the vicinity of 
xi given by 

. 1 d' 3>i 

a = ai — i- — 

dx 

The real part of this effective wave number represents the phase change of the 
disturbance while the imaginary part depicts the growth rate, both corresponding 
to the quantity chosen. A disturbance (’Fj) is unstable if the imaginary part 
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is less than zero. The updating procedure of a is repeated by using (17) until the 
change in a is smaller than a prescribed tolerance (typically 10 -12 ). 

Since the shape function vector is a function of y and contains five depen- 
dent variables (/3,u, etc.), the updating procedure above is equivalent to choosing 
a normalization of the disturbance vector such that d'&i/dx is zero at a particular 
y location. Accordingly, the value of o computed by (17) will depend on the y 
coordinate and the selected dependent variable \&i. In this study, we have used the 
shape function u (or t for compressible flows) at various y locations or the energy 
integral ( E = f^°p(u 2 + v 2 + w 2 )dy), which is independent of the y coordinate, 
to update the wave number a and the resulting non-parallel growth rate (which 
also depends on the dependent variable and y coordinate chosen to measure the 
growth rate) appears to be very weakly dependent upon the normalization chosen 
(see discussion in section IV). 

The solution of (14) requires proper boundary conditions in the wall-normal 
direction. We apply the homogeneous Dirichlet conditions 

u=v = w = T = 0, y = 0 (18) 

at the wall and in the free-stream 

u — v = w = T — 0, y — » oo; (19) 

although, these can be easily replaced by other conditions such as the Rankine- 
Hugoniot conditions at the shock (Chang et al., 1990) for supersonic flows. Non- 
homogeneous boundary conditions can also be imposed. 

Nonlinear PSE 


In the linear PSE approach described above, the disturbance amplitude is as- 
sumed to be infinitesimally small so that the nonlinear interaction of waves with dif- 
ferent frequencies and spanwise wave numbers is neglected. When finite- amplitude 
waves are present in the flow, the linear approach is no longer valid. For nonlinear 
studies, we assume that the total disturbance is again periodic in time and in the 
spanwise direction, thus, the total disturbance function <f> can be expressed by the 
following Fourier series 


OO OO r x 

, \ — ' _ , . i( / a mn (x)dx+nPz-mwt) 

4> = 2^ *mn(z,y)e 


771= —oc n= — OO 


( 20 ) 


where a mn and are the Fourier components of the streamwise wave number 

and shape function corresponding to the Fourier mode (mu,n0). The frequency 
uj and wave number /? are chosen such that the longest period and wave length 
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are 2v/u and 27r//3 in the temporal and spanwise domains, respectively. For most 
stability problems of interest, it is sufficient to truncate (20) to only a finite number 
of modes 


M N f . 


m= — M n~ — N 


( 21 ) 


where M and N are the total number of modes kept in the truncated Fourier series. 
For all nonlinear results presented in this study, we apply both the temporal and 
spanwise symmetry conditions whenever applicable, i.e., only one quarter of modes 
(m ranging from 0 to M and n ranging from 0 to N) are computed in the marching 
process. 

We now substitute Eq. (21) into our nonlinear governing equation (6) and 
perform harmonic balance (collect terms with the same spanwise wave number and 
frequency) for both linear and nonlinear terms. The resulting governing equations 
for the shape function of a single Fourier mode (m,n) become 


t) +A ^ mn . £ d%n. 

^mn * mn r rv i -°mn r\ 

ox oy 


v i d 2 * mn 
v vv Q y 2 


+ F mn JA 


mn 


(22) 


where matrices D mn , A mn and B mn are given by 

D m n = -imu>F l + D 1 + ia mn A l + inftC 1 

~ ” a mn) v L + na mn /3V‘ z 

+ n 2 /3 2 V zz 

Amn — A 2 ioL mn V xx in(3V xz 

B m n — B i&mnV X y 


and the quantity A mn is 


i f* a rnn (x)d2 


The nonlinear forcing function F mn is the Fourier component of the toteil forcing 
defined by Eq. (7) and can be evaluated by the Fourier series expansion of F n 


M N 

F"(i, »,.,<)= E E (23) 

m= — jVf n= — JV 


The Fourier decomposition of Eq. (23) can be done by using the Fast Fourier 
Transform (FFT) of F n , which is evaluated numerically in the physical space. In 
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equation (22), a parabolizing procedure similar to that used in the linear PSE has 
been employed in order to obtain a marching solution. 

As in the linear PSE, the determination of the wave number a mn plays an 
important role in maintaining numerical stability. The procedure described above 
for computing a for linear disturbances can also be used for the determination of 
a mn . However, when all Fourier modes are nearly phase-locked (as is evident when 
parametric resonance of secondary instability takes place, see e.g. Kachanov and 
Levchenko (1984)), one may assume that the wave number is given as 


®mn — (ma r , <T mn ) 


where a r is the real part of Qio and cr mn denotes the growth rate of the mode (m, n). 
Each mode can have a different imaginary part <r mn , while the real part is updated 
according to the phase change of the dominant fundamental mode. Additional 
phase shifts in the harmonics are included in the evolution of the shape functions of 
the harmonic waves; therefore, all Fourier modes are not necessarily phase-locked. 
It needs to be pointed out that the “nearly” phase-locking assumption (since the 
evolution of shape functions may shift the phase slightly) mentioned above is used 
for convenience and is not necessary for the nonlinear analysis using PSE. In a 
later section, we will provide an example of a nonlinear calculation where we let 
the disturbances evolve with and without the phase- locking rule. Use of the phase- 
locking rule, when applicable, saves computational cost. 

The nonlinear PSE for a single Fourier mode, equation (22), is equivalent to 
the linear PSE given in (14) with a frequency mw and a spanwise wave number nf3 
with the addition of a forcing function. Since the forcing function acts as a “source 
term” of the equation, the boundary conditions and solution procedure described 
above for the linear PSE can be directly applied to the nonlinear system, except for 
the modes with zero frequency (m = 0). These zero frequency modes are denoted 
as the mean flow correction (if n = 0) or longitudinal vortex modes (if n ^ 0). For 
these modes, as in Gortler vortex problem, the pressure gradient dp/dx drops out 
making the equations fully parabolic. 

The boundary conditions given in Eqs. (18) and (19) can be applied to the 
longitudinal vortex mode without modification. For the mean flow correction, the 
free-stream conditions are replaced by 


^oo 


dVQQ 

dy 


woo — Too = 0, 


y — * oo 


(24) 


to account for the change of displacement thickness due to the correction of the 
mean flow profile (u + tioo) arising from nonlinear interactions. This Neumann 
condition for the normal velocity allows the mean flow given by the boundary-layer 
solution to adjust itself in order to assure mass balance. 


11 


Ill Numerical Procedure 


In this paper, we only consider the compressible stability of two-dimensional 
boundary-layer flow past a flat plate. The mean flow solution is obtained by solv- 
ing the self- similar boundary layer equations. By using the Mangier- Levy- Lees 
transformation, the boundary-layer equations are transformed into a set of ordi- 
nary differential equations (ODE’s). A fourth-order compact scheme is employed to 
solve these ODE’s. Details of the numerical procedures are given in Malik (1990b) 
and will not be repeated here. 

Numerical solution of the parabolized stability equations (14) or (22) requires 
discretization in both x and y directions. Since the boundary layer grows in the 
streamwise direction, we expect that the solution for the shape functions will also 
grow. To ensure sufficient resolution as the disturbances evolve downstream, dis- 
cretization in the wall-normal direction must be able to account for the growth of 
the boundary layer. Instead of solving equations (14) and ( 22 ) in Cartesian coor- 
dinates, we transform these equations to a generalized coordinate system defined 

by 

v = vi x ,y) 

in order to facilitate numerical computations on a “growing mesh” or curved wall 
geometries. After this transformation, Equation (14) becomes 



~dV ~ d 2 V 

D * + A a, I +B ej = v "W 


where the coefficient matrices are given by 


(26) 


D = D 

A = ( X A + i y B 

V 1 d r ? 2 

B = v ,A + „B--f^ £) 

Vrjrf — VyVyy 

The Jacobian of the transformation J is defined as 


J — ix'Hy 

Equation (22) can be transformed in a similar fashion. 

Using transformation (25), we map the computational grid into a uniform mesh 
with constant increments in £ and 77 coordinates. For most of our calculations, we 
use a constant step size in x while the grid is clustered near the wall to resolve the 
rapid change inside the boundary layer. For high Mach number calculations, we also 
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cluster the grid near the critical layer located near the boundary-layer ed ge. The 
stretching along the y direction is based on the local length scale l x (l x = \Jv t xju t ) 
which increases with the boundary-layer growth. The same grid distribution based 
on y/l x is used for all x locations while l x increases with x. 

We use indices i and j to denote the grid index along the streamwise (x) and 
wall-normal (y) directions, respectively. In the streamwise direction, we use the 
second-order backward difference 


d* 


(3*i j - + *i_ w )/2A{ 


for all x locations, 
ence is employed, 
is then 


except for the starting plane where a first-order backward differ- 
The resulting discretized equation for the i-th streamwise plane 


[t>+ 2 

i-i j - 


(27) 


In the wall-normal direction, we employ a fourth-order accurate finite-difference 
scheme. The two-point fourth-order scheme by Malik et al. (1982) is also used for 
normal derivatives; however, it requires that the normal mean velocity v is non- zero 
and hence is not generally applicable for all problems. 

For the uniform mesh in the £ — r) plane, the normal derivatives in Eq. (26) 
are discretized according to the following fourth-order central difference formulae: 


drj 12A»7 


~ q ^2 =(-^«J+2 + 16^«, >+i - 30tfj,j 
+ - V / 12 Ar} 2 . 

For the grid point next to the boundary, the above five-point scheme is replaced 
with the second-order scheme 


drj 




^ - 2 + * U -i)/2W- 

At the boundary, five boundary conditions are needed for five dependent variables 
in The no-slip and free-stream boundary conditions given in (18) and (19) 
are used. In addition, we apply the discretized continuity equation as the fifth 
boundary condition both at the wall and the free-stream. Substituting the above 
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normal derivatives into (27) for all interior points and coupling with the boundary 
conditions result in a block penta-diagonal system of equations at each x location 
with a block size of 5 x 5. This block matrix can be solved by the standard LU 
decomposition method. 


IV. Results and Discussion 

To demonstrate the capability of the PSE approach, we perform both linear 
and nonlinear calculations for various Mach numbers. In the linear results, the 
main focus will be on the non-parallel effect, and in the nonlinear regime, PSE 
calculations are carried all the way to the early stage of transition. 

In the following discussion, we define the growth rate in non-parallel boundary 
layers according to Eq. (17), i.e., for any given flow variable ip ( for instance, p, u, 
etc.), the growth rate a is defined as 

a = -Im(a) + Re(^). (28) 

The second term on the right hand side of the above equation is a function of y\ 
therefore, the growth rate in a non-parallel boundary layer depends upon the dis- 
tance normal to the wall. We note here that although Eqs. (17) and (28) are derived 
based on the same concept, they have different physical meanings. Briefly, Eq. (17) 
is used to normalize the disturbance vector and determine the wave number a. For 
each normalization, corresponding to different Ski chosen in (17), the growth rate 
for any given variable at any y location is to be evaluated using Eq. (28). For the 
results presented herein, we compute the growth rate at the corresponding loca- 
tion where the fluctuation reaches its maximum value or based on the disturbance 
kinetic energy integral, 

<te = —Im(ot) + —(lny/E) 

where E is defined by 

1*00 

E = I ( u 2 + v 2 + w 2 )dy 

Jo 

for the incompressible limit and by 

J v*oo 

' p(u 2 + v 2 + w 2 )dy 
o 

for general compressible flows. In supersonic wind tunnel experiments, the growth 
rate is usually measured for the mass flow fluctuation. We define the mass flow 
fluctuation as 

(pu)' = pu + pu 1 . (29) 
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Linear PSE 


As mentioned in the previous section, the streamwise wave number a depends 
upon the variable 'k i chosen and the y location where (17) is applied. To demon- 
strate that the resulting non-parallel growth rate is very weakly dependent upon 
various normalizations, we first perform calculations for a Mach 1.6 boundary layer 
by using different dependent variables to update a. These variables include u, T 
(evaluated at various y locations as shown in the figure) and the kinetic energy inte- 
gral E defined above. Figure 1 shows the resulting imaginary part of the converged 
value of a for the various norms. The results reveal that a strongly depends on 
the norm chosen. The corresponding effective growth rates cr pu , ax and a e, eval- 
uated by using (28) at their maximum locations (for a pu and ax only) are shown 
in Fig. 2. It shows that the total growth rates depend on how they are measured 
(for instance, ax and oe are different); however, each non-parallel growth rate (e.g. 
a pu ) appears to be independent of the normalization procedure because results from 
various norms collapse into one single curve. Similarly, although not shown here, 
the non-parallel wave number, evaluated by 


a 


np 


_ . . ,1 dip . 

= Re(a) - 


is also weakly dependent on the normalization. The above results indicate that 
although different norms result in different values of a, the total growth rate (and 
wave number) by accounting for the evolution of shape function in the streamwise 
direction remains the same regardless of the norms. For the results presented herein, 
we use the kinetic energy integral E in (17) to update a. 

To verify the numerical algorithm, the first test case studied is an incom- 
pressible flow case. The incompressible results were obtained by choosing a Mach 
number of 10 -6 in our compressible formulation. Linear non-parallel results are 
available for incompressible boundary layer flow by using local methods from many 
authors(e.g., Gaster, 1974). The neutral points obtained from our PSE calcula- 
tions agree very well with those from Gaster’s (1974) non-parallel method. Figure 
3 shows the computed variation of the growth rates (<r„, a v and <je ) with Reynolds 
number ( R = y/u^xju^) for a representative non-dimensional frequency ( F = uiR) 
of 1.12 x 10 -4 . The growth rates from multiple-scales method are shown by symbols. 
The results shown in Figure 3 reveal that the neutral curve near the upper branch is 
shifted to higher Reynolds numbers due to non-parallel effect as was found by Gaster 
(1974) and linear PSE results agree very well with those from the multiple-scales 
approach. 

The second test case was chosen to be the Mach 1.6 case studied by El-Hady 
(1991) using the multiple- scales approach. The frequency was fixed at 0.4 x 10~ 4 
and variable transport properties were used. Calculations were performed for both 
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2-D and 3-D linear disturbances with an oblique wave angle of about 50° for the 
latter. The growth rate of the mass flow fluctuations (defined in Eqs. (28) and 
(29)) from our PSE calculations together with the multiple-scales results are plotted 
along with the growth rates obtained by quasi-parallel linear stability theory in Fig. 
4. Our PSE results agree quite well with those obtained from the multiple-scales 
approach. The results also indicate that for the first mode disturbance at Mach 
1.6, flow non-parallelism has more effect on three-dimensional disturbances than 
on two-dimensional ones. Results obtained at higher Mach numbers also show a 
noticeable non-parallel effect on the first-mode instability. 

We now show some results for the Mach 4.5 flat-plate flow, which is subject 
to second-mode instability (Mack 1984). Calculations were performed for a dis- 
turbance frequency of F — 1.2 x 10 -4 with different streamwise resolutions. We 
used step sizes ranging anywhere from 64 steps per wavelength to only one step 
per wavelength. The results for the second-mode growth rate based upon the to- 
tal kinetic energy are plotted in Figure 5. There is essentially no difference in the 
growth rate results when two or more steps per wavelength are used. The reason 
why only two points per wavelength could yield such accurate growth rates lies in 
the fact that most of the wave information is absorbed in the complex wavenumber 
a. In contrast, direct numerical simulation (DNS) of Navier-Stokes equations would 
require many more points per wavelength for comparable accuracy. 

To further verify our linear results, we compare non-parallel evolution of a 
second mode disturbance with a frequency of 2.2 x 10 -4 with DNS. In Fig. 6, 
the maximum amplitudes of various flow quantities are plotted against Reynolds 
numbers for both PSE and DNS. The PSE results obtained by using only 7 steps 
per wave length agree very well with DNS results using 16 steps per wave length. 
The PSE calculation took about 100 seconds CPU time while the DNS required 
more than 40 hours on a Cray-YMP. Details of the comparison including nonlinear 
disturbances and the spatial DNS algorithm are given in Pruett and Chang (1993). 

Nonlinear PSE 


a. Computation of o 

In a previous section we mentioned that the wavenumber a for the harmonics 
may be determined either by the phase-locking rule or by using Eq. (17). It is known 
that the nonlinear wave interaction is dependent on the phase-difference between 
various modes. Therefore, it is essential that nonlinear PSE approach must not 
require phase-locking as a fundamental assumption; although, it may be used as a 
convenience for problems where phase-locking happens anyway. 

In order to demonstrate that phase-locking is not a basic assumption for non- 
linear PSE computations, the following test has been performed. Nonlinear calcula- 
tions have been done for a flat-plate boundary layer in the incompressible limit. A 
two-dimensional wave with a frequency F — .86 x 10 -4 and an initial amplitude of 
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.25% at R = 400 is introduced in the boundary layer and the evolution of this wave 
along with its various harmonics is monitored. Calculations were performed in two 
different ways. First, the wavenumber a was computed for the fundamental wave 
according to Eq. (17) and the phase-locking rule was used for all the harmonics. In 
the second set of calculations, wavenumbers for the fundamental and all the harmon- 
ics were computed independently by using Eq. (17). The computed results for the 
amplitude of u velocity (fundamental, its four harmonics and meanflow correction) 
are presented in Fig. 7(a). It can be seen that only very minor differences appear 
between the two sets of calculations and these also tend to disappear as the calcu- 
lations are marched away from the inflow boundary. The second set of calculations 
takes about 50% more computer time due to the iterations involved in determining 
a for the wave harmonics. Hence, it is expedient to use the phase-locking rule for 
problems where this may be the outcome in any case. 

In Fig. 7(b), the same nonlinear PSE results are compared with the spatial 
incompressible DNS results for fundamental, first harmonic and the mean flow dis- 
tortion modes. Both PSE and DNS start with the same initial conditions, i.e., a fun- 
damental disturbance (1,0) at R = 400 and all harmonic waves including the mean 
flow distortion are generated through nonlinear interactions. The good agreement 
between DNS and nonlinear PSE indicates that the parabolizing approximation in 
the PSE approach does not introduce any severe error and all detailed nonlinear 
features are properly captured. Details of the comparison including disturbance 
profiles can be found in Joslin et al. (1992). 
b. Second-Mode Instability at Mach 4.5 

To verify the nonlinear PSE algorithm, we choose the nonlinear second mode 
simulation at Mach 4.5 investigated by Erlebacher and Hussaini (1990) using the 
temporal DNS approach. As in the temporal DNS approach, we assume that the 
mean flow is parallel and study the spatial evolution of disturbances in the presence 
of nonlinear interactions. The initial conditions were provided by the eigensolution 
from the linear theory at R = 781 and four Fourier modes ( M = 3) were kept in the 
truncated series. In our PSE calculation, the disturbance is assumed to be periodic 
in time and the nonlinear evolution is carried downstream in x as opposed to the 
temporal DNS approach where the disturbance is periodic in x and integration is 
carried in time. It was found in Erlebacher and Hussaini (1990) that due to non- 
linear effect, the growth rate of the fundamental disturbance strongly depends on 
y and there exists a sharp decrease in the local growth rate near the critical layer. 
The growth rates based on t<io from our PSE results are shown in Figure 8(a) for 
different x locations. The growth rate is initially uniform at the starting location 
(x = 0A). As the disturbances are evolving downstream, nonlinear effects observed 
in Erlebacher and Hussaini (1990) are evident in the present spatial calculations. 
For comparison, their temporal DNS results are shown in Figure 8(b) at different 
time levels represented as multiples of the temporal period r. Figure 9 depicts the 
amplitudes of the density fluctuation of the first harmonic for both PSE calcula- 
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tions and the DNS results. The DNS results in Figure 9 are re-scaled to facilitate 
comparison. Qualitatively, all nonlinear features observed in DNS, including the 
kink near the boundary-layer edge, are properly resolved in our PSE results, 
c. Subharmonic and Fundamental Resonance at Mach 1.6 

Numerical simulation of incompressible flows have shown that the rapid growth 
of three-dimensional secondary disturbances is followed by breakdown to turbulence. 
Secondary instability is triggered when the primary disturbances reach sufficiently 
high amplitudes. To show the capability of the PSE approach in simulating transi- 
tion onset, we also perform nonlinear calculations to study the secondary instability 
mechanism. We carry our calculations all the way to transition for both K-type (fun- 
damental) and H-type (subharmonic) breakdown. We choose Mach 1.6 flat plate 
flow with a primary disturbance frequency of 0.5 x 10 -4 . The same flow conditions 
were also used by Thumm et al. (1989) in their spatial Navier-Stokes simulations 
of a compressible boundary layer. 

We first perform a series of calculations to determine the amplitude of the 
primary disturbance which will trigger the secondary instability. The PSE calcu- 
lation is initiated at a Reynolds number R of 460 where we impose a primary 2-D 
wave (mode (2, 0)) obtained by a local eigenvalue calculation and two subharmonic 
waves ((1,1) and (1,-1) modes) by using the compressible secondary instability 
theory (Ng and Erlebacher, 1992) and all the remaining harmonics are assumed to 
have zero amplitudes. Initial amplitudes of the primary disturbances are set to be 
3%, 1.1% and 0.6% at the inflow plane which corresponds to the 5%, 2% and 1% 
(the maximum amplitudes near the vibrating ribbon) cases given in Thumm et al. 
(1989). The initial amplitudes of the subharmonic waves are fixed at 0.019% for 
all three cases. The spanwise wave number of the subharmonic mode is fixed at 
fi/R = 0.53 x 10 -4 which corresponds to an oblique wave angle of 45°. Six temporal 
Fourier modes and three spanwise modes (M = 5 and N = 2) are kept in the Fourier 
series. The evolution of both primary and subharmonic disturbances are shown in 
Figure 10. Qualitatively, our results agree with those of Thumm et al.(1989). Any 
quantitative differences are due to different initial conditions. We find that a 1.1% 
initial amplitude (2% case in Thumm et al. (1989)) for the primary mode is enough 
to trigger the secondary growth; however, the onset of secondary growth for this 
case occurs at R = 800 where the primary wave is about to decay. We continue 
the PSE calculations beyond R = 1050 for this case and find that the secondary 
disturbance eventually saturates and the flow does not reach the transitional stage. 

To carry the 3% case to the transition stage, we made another calculation with 
more Fourier modes (M = 7 and N = 4) and the maximum velocity amplitudes of 
some representative modes are given in Figure 11. Besides the fundamental mode 
(2,0) and the subharmonic mode (1,1), higher harmonics in time and spanwise 
domain are also excited due to nonlinear interaction. Initially, the (4, 0) mode gains 
energy from self-interaction of the (2, 0) mode and the interaction of (2, 0) and (1, 1) 
produces the (3, 1) mode. When the subharmonic mode grows due to the onset of 
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secondary instability, its harmonic (2,2) also grows at slightly higher rate. The 
streamwise vortex mode (0,2) arises due to the interaction of (1, 1) mode and its 
complex conjugate ( — 1,1). As all these modes continue to grow, more and more 
modes are excited. The energy cascade exhibits a staggered pattern. For instance, 
among the two-dimensional modes, only (2,0), (4,0), (6,0), etc. gain energy; while 
for 1/3 modes, only (1, 1), (3, 1), (5, 1), etc. are excited. The remaining modes (e.g. 
(1,0), (3,0), (0,1), (2,1) etc.) remain unexcited throughout the calculation. The 
above staggered energy cascade is typical for subharmonic secondary instability. 
The secondary amplitude overtakes the primary at about R = 980 and reaches an 
equilibrium state around R = 1100. At this stage, many harmonic waves reach 
fairly high amplitudes as the flow heads for transition. We plot the time sequence 
of spanwise vorticity contours at the peak and valley planes (corresponding to the 
maximum and minimum disturbance rms amplitudes, respectively) in Figures 12(a) 
and 12(b). As can be seen, the vorticity pattern doubles its wavelength for x > 
2200 (x is normalized w.r.t. the boundary layer length scale l at the initial plane) 
indicating the presence of high-amplitude subharmonic wave. It is also evident 
that the vortex roll-up results in a distinct kink in the shear layer. Towards the 
end of the computation, regions of intense vorticity near the wall begin to appear 
indicating that flow is heading for breakdown. Figure 13 shows the streamwise 
velocity contours in the x-z plane for a wall normal distance of y = 2.3, where the 
TS wave reaches its maximum according to the linear solution. The flow is initially 
two-dimensional and three-dimensional effect becomes important for x > 1800. For 
x > 2000, a staggered contour pattern is evident. This pattern is associated with 
the lambda vortex structure, a distinct characteristic of subharmonic breakdown, 
as observed in many incompressible experiments (e.g. Corke and Mangano (1989)). 

Nonlinear PSE calculations are also performed for the same Mach 1.6 case but 
for a fundamental- type secondary resonance. The initial amplitude of the primary 
wave is again 3% and that of the secondary is taken to be 0.005% to minimize 
nonlinear interaction close to the starting location. The spanwise wave number is 
P/R = 1.52 x 10 4 (oblique wave angle of 60° for the secondary wave) and the pri- 
mary wave frequency is again 0.5 x 10 -4 . The initial conditions for our marching cal- 
culation consist of a 2-D primary wave (mode (1, 0)), two oblique fundamental-type 
secondary disturbances (mode (1,1), (1,-1)) and the longitudinal vortex (mode 
(0, 1)). The same number of Fourier modes as in the subharmonic case is used. 

Nonlinear evolution of the maximum rms amplitude of u' is shown in Figure 14. 
Initially, the dominant modes are (1,0), (1, 1), (0, 1) and (2, 0) (the first harmonic of 
the fundamental 2D mode). Unlike the subharmonic case, all harmonic waves (both 
odd and even modes) gain energy directly from nonlinear interaction. Among them, 
the (2,1) (due to (1,0) and (1,1)) and (1,2) (due to (0,1) and (1,1)) modes are 
more noticeable. For Reynolds numbers beyond 870, the. spectrum is rapidly filled 
with high-amplitude disturbances and the flow is heading towards transition. As 
compared to the subharmonic breakdown, transition location shifts upstream due 
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to the larger growth rate of the secondary disturbance as a consequence of higher 
oblique wave angle. 

The time sequence of spanwise vorticity contours over a period of the primary 
wave is shown in Figures 15(a) and 15(b) for the peak and valley planes, respec- 
tively. In contrast to the subharmonic breakdoown, the wave length remains the 
same throughout the whole computational domain. One important characteristic of 
the K-type breakdown is the appearance of aligned lambda vortices. This is better 
visualized in the streamwise velocity contours shown in Figure 16 for x > 1300. 
Similar to that observed in incompressible simulations of Zang and Krist (1989), 
regions of intense shear begin to appear near the end of the computational domain 
in Figures 15(a) and 15(b). This indicates that flow has just entered the transitional 
stage. It is confirmed by plotting the average wall shear in Figure 17. The com- 
puted wall shear is slightly above the laminar value for most of the computational 
domain. Only towards the end, wall shear significantly departs from the laminar 
value indicating the onset of transition. In this way, PSE provides the prediction 
of boundary-layer transition for the imposed initial conditions. The PSE wall shear 
lies above the laminar value right from the beginning because of the relatively high 
amplitude of the 2-D primary disturbance needed for transition in supersonic flow. 
Since most amplified waves in supersonic flow are not two-dimensional, oblique 
primary modes may lead to transition for lower initial amplitudes. In order to 
carry the calculations further into transitional regime, more spanwise and temporal 
resolution will be required. It remains to be seen how far PSE can proceed into 
the transitional zone. The computational time used for the results presented in 
Figures 14-17 was 15 minutes on a Cray-YMP machine. Similar results from full 
compressible Navier- Stokes equations would require 0(50) hours. 


V. Conclusions 

Linear and nonlinear compressible boundary-layer stability is studied by using 
the PSE approach. Several issues concerning the characteristics of the paraboliza- 
tion and the updating of the streamwise wave number are also discussed. The 
governing equations are solved by using second-order backward differences for the 
streamwise derivatives while the wall-normal direction is discretized by a fourth- 
order accurate finite- difference scheme. 

Non-parallel flow effects have been studied for linear disturbances. For oblique 
waves of the first mode type, the departure from the parallel results is more pro- 
nounced as compared to that for the two-dimensional waves. Our linear results are 
in good agreement with those from the multiple-scales approach, as well as those 
from full Navier- Stokes equations. 

Nonlinear PSE calculations are carried all the way to the early stage of tran- 
sition for a Mach 1.6 flow. Both the subharmonic and fundamental types of break- 
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down are studied by the current PSE approach. Qualitatively, these breakdown 
processes are similar to the ones in incompressible boundary layers, except that 
high amplitudes of the 2-D primary wave are required. The promising results of 
our PSE calculations show that this new approach is a powerful tool for the study 
of boundary-layer stability and transition prediction. The parabolized form of the 
governing equations allow the numerical solution to be obtained in a computational 
time which is orders of magnitude lower than that required for direct simulation of 
Navier- Stokes equations. 
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Figure 1. converged values of — a; by using various normalizations for a Mach 1.6 
flow with F = 0.5 x 10~ 4 
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Figure 2. Total non-parallel growth rates <7 pu , a T and a E for different normaliza- 
tions corresponding to Figure 1. 
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Figure 5. Growth rates based on kinetic energy for a Mach 4.5 second-mode distur- 
bance with F — 1.2 x 10 -4 using various streamwise step sizes. 



Mach 4.5 F=2.2X10 



Figure 6. Comparison of linear PSE and spatial DNS results for a Mach 4.5 second- 
mode disturbance with F = 2.2 x 10~ 4 



incompressible nonlinear F=86 



Figure 7. Nonlinear evolution of a fund am ental 2-D wave with F = 0.86 x 10 -4 in 
the incompressible limit: (a)showing the velocity amplitudes both with and without 
phase- locking rule (b)compaxing with spatial DNS results. 
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Figure 8. Nonlinear second mode simulation on a parallel mean flow of Mach 4.5 
and R = 781, showing the evolution of growth rates against y at various downstream 
distances represented as multiples of the fundamental wave length A (for PSE) or 
the temporal period r (for DNS). (a)Nonlinear PSE (b)DNS 
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Figure 9. Nonlinear second mode simulation on a parallel mean flow of Mach 4.5 and 
R - 781: showing the evolution of density profiles against y for the first harmonic 
wave (mode (2,0) ). Lines are from nonlinear PSE calculations and symbols are 
from DNS. 





Mach 1 .6 Subharmonic Breakdown F=50 








Figure 12. Time sequence of spanwise vorticity contours for the Mach 1.6 subhar- 
monic breakdown (x and y are streamwise and wall-normal coordinates normalized 
by the boundary-layer length scale at R = 460): (a) Peak plane, (b) Valley plane. 
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Fig. 12(b) 
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Mach 1 .6 Subharmonic Breakdown 
U Velocity Contours at y/l 0 =2.3 



Figure 13. Streamwise velocity contours in x — z plane at y = 2.3 for the subhar 
monic breakdown. 








Mach 1 .6 Fundamental Breakdown F=50 
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Figure 14. Amplitude evolution of various Fourier modes for a Mach 1.6 fundamen- 
tal breakdown (F = 0.5 x 10 -4 and ft/R = 1.52 x 10 -4 ) with initial amplitudes of 
3% and 0.005% for the primary and secondary disturbances at R = 460. 




Figure 15. Time sequence of spanwise vorticity contours for the Mach 1.6 funda- 
mental breakdown : (a) Peak plane, (b) Valley plane. 
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Fig. 15(b) 
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Mach 1 .6 Fundamental Breakdown 
U Velocity Contours at y/l x =2.3 
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Figure 16. Streamwise velocity contours in x — z plane at y — 2.3 for the fundamental 
breakdown. 
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